rm(list=ls())

load("TSCS_data.RData")

set.seed(0211)

models_plm <- list()

for (iteration in 1:100) {
  
  panel %>% mutate(antipluralism_gov_seat_share_lag_draw = effectsize::standardize(lag(panel$antipluralism_gov_seat_share,sample(1:2,1)),two_sd = T),
                   citizen_support_lag_draw = effectsize::standardize(lag(panel$citizen_support,sample(1:2,1)),two_sd = T)) ->
    panel
  
  model <-   plm(
    downturn ~ lag(downturn,1:2)  +
      antipluralism_gov_seat_share_lag_draw*citizen_support_lag_draw + 
      lag(polarization,sample(1:2,1)) +
      lag(clientelism_gov_seat_share,sample(1:2,1)) + 
      lag(region_citizen_support,sample(1:2,1)) + 
      lag(region_antipluralism_gov_seat_share,sample(1:2,1))+ 
      lag(region_liberal_democracy,sample(1:2,1)) + 
      lag(democratic_stock,sample(1:2,1)) + 
      lag(GDP_capita,sample(1:2,1)) +
      lag(life_expectancy,sample(1:2,1)),
    data = panel,
    model = "within",
    effect = "twoways")
  
  models_plm[[iteration]] <- model
}

coeffs_lm <- vector()
se_lm <- vector()

for (i in 1:100) {
  coeffs_lm[i] <- models_plm[[i]][["coefficients"]][["antipluralism_gov_seat_share_lag_draw:citizen_support_lag_draw"]]
}

vcm_lm <- lapply(models_plm, function(x) plm::vcovBK(x, cluster=c("group")))
rse_lm <- lapply(vcm_lm, function(x) sqrt(diag(x)))

for (i in 1:100) {
  se_lm[i] <- rse_lm[[i]][["antipluralism_gov_seat_share_lag_draw:citizen_support_lag_draw"]]
}

coeff_df <- as.data.frame(cbind(coeffs_lm,se_lm))
coeff_df %>% mutate(margin_of_error = se_lm*1.63,
                    CI_lower = coeffs_lm - margin_of_error,
                    CI_upper = coeffs_lm + margin_of_error) %>%
  arrange(coeffs_lm)  %>%
  mutate(model = row_number())->
  coeff_df

pdf(file = "fig_e2a.pdf",width = 3,height = 9)

ggplot(coeff_df,aes(x=coeffs_lm,y=model))  + 
  geom_errorbarh(aes(xmin=CI_lower, xmax=CI_upper), colour="grey") +
  theme() + xlab("Coefficient") + ylab("Model")  +
  geom_vline(xintercept = 0) + geom_point() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 

dev.off()

# GMM

set.seed(0211)

models_pgmm <- list()

for (iteration in 1:100) {
  
  panel %>% mutate(antipluralism_gov_seat_share_lag_draw = effectsize::standardize(lag(panel$antipluralism_gov_seat_share,sample(1:2,1)),two_sd = T),
                   citizen_support_lag_draw = effectsize::standardize(lag(panel$citizen_support,sample(1:2,1)),two_sd = T)) ->
    panel
  
  model <- pgmm(
    downturn ~ lag(downturn,1:2)  +
      antipluralism_gov_seat_share_lag_draw*citizen_support_lag_draw + 
      lag(polarization,sample(1:2,1)) +
      lag(clientelism_gov_seat_share,sample(1:2,1)) + 
      lag(region_citizen_support,sample(1:2,1)) + 
      lag(region_antipluralism_gov_seat_share,sample(1:2,1))+ 
      lag(region_liberal_democracy,sample(1:2,1)) + 
      lag(democratic_stock,sample(1:2,1)) + 
      lag(GDP_capita,sample(1:2,1)) +
      lag(life_expectancy,sample(1:2,1)) |
      lag(downturn, 3:5),
    panel,
    effect = "individual",
    model = "onestep",
    transformation = "ld",
    indexes = c("country", "year"),
    robust = T)

  models_pgmm[[iteration]] <- model
}

coeffs_gmm <- vector()
se_gmm <- vector()

for (i in 1:100) {
  coeffs_gmm[i] <- models_pgmm[[i]][["coefficients"]][["antipluralism_gov_seat_share_lag_draw:citizen_support_lag_draw"]]
}

vcm_gmm <- lapply(models_pgmm, function(x) plm::vcovHC(x, cluster=c("group")))
rse_gmm <- lapply(vcm_gmm, function(x) sqrt(diag(x)))

for (i in 1:100) {
  se_gmm[i] <- rse_gmm[[i]][["antipluralism_gov_seat_share_lag_draw:citizen_support_lag_draw"]]
}

coeff_df <- as.data.frame(cbind(coeffs_gmm,se_gmm))
coeff_df %>% mutate(margin_of_error = se_gmm*1.63,
                    CI_lower = coeffs_gmm - margin_of_error,
                    CI_upper = coeffs_gmm + margin_of_error) %>%
  arrange(coeffs_gmm)  %>%
  mutate(model = row_number())->
  coeff_df

pdf(file = "fig_e2b.pdf",width = 3,height = 9)

ggplot(coeff_df,aes(x=coeffs_gmm,y=model))  + 
  geom_errorbarh(aes(xmin=CI_lower, xmax=CI_upper), colour="grey") +
  xlab("Coefficient") + ylab("Model")  +
  geom_vline(xintercept = 0) + 
  geom_point() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 

dev.off()